The effect of body weight on the knee joint biomechanics based on subject-specific finite element-musculoskeletal approach

Knee osteoarthritis (OA) and obesity are major public health concerns that are closely intertwined. This intimate relationship was documented by considering obesity as the most significant preventable risk factor associated with knee OA. To date, however, the effects of obesity on the knee joint's passive-active structure and cartilage loading have been inconclusive. Hence, this study investigates the intricate relationship between obesity and knee OA, centering on the biomechanical changes in knee joint active and passive reactions during the stance phase of gait. Using a subject-specific musculoskeletal and finite element approach, muscle forces, ligament stresses, and articular cartilage contact stresses were analyzed among 60 individuals with different body mass indices (BMI) classified under healthy weight, overweight, and obese categories. Our predicted results showed that obesity significantly influenced knee joint mechanical reaction, increasing muscle activations, ligament loading, and articular cartilage contact stresses, particularly during key instances of the gait cycle—first and second peak loading instances. The study underscores the critical role of excessive body weight in exacerbating knee joint stress distribution and cartilage damage. Hence, the insights gained provide a valuable biomechanical perspective on the interaction between body weight and knee joint health, offering a clinical utility in assessing the risks associated with obesity and knee OA.


Power analysis
Using data ranges previously reported in the literature for outcomes treating the impact of obesity on the biomechanics of the lower extremities during gait [1-5], a preliminary power analysis was performed using G*Power [6,7] which gave a conservative estimate varying from 6 to 18 subjects in each group to be sufficient to detect meaningful differences between groups with a power of 1 -β = 0.80.This estimate is consistent with sample sizes reported for similar studies [8,9].

Inclusion Criteria
The criteria for participant inclusion were as follows: only males between the ages of 18 and 40 years were eligible.To qualify, participants must have maintained a stable body weight, with a variation of less than 2.5 kg in the last three months.They also need to engage in less than 30 minutes of moderate physical activity each week and not more than three days per week.
Additionally, candidates with no previous history of knee pain or lower limb surgeries were considered.Finally, only those who gave informed consent were included in the study.

Exclusion Criteria
Exclusion criteria for the study were as follows: participants who had a history of knee pain or had undergone any lower limb surgeries were not eligible.Also, those who participated in moderate physical activity for over 30 minutes on at least three days per week were excluded.Individuals whose body weight had shifted by more than 2.5 kg within the previous three months were also ruled out.Additionally, participants displaying signs of osteoarthritis in their knee joints, as indicated by X-ray images, were excluded from the study.

Subject population
Based on the power analysis, inclusion, and exclusion criteria outlined previously, the research recruited 60 male participants who were similar in age, daily activity levels, and body mass stability (as shown in Table 1 of the manuscripts).The participants were recruited via announcements and advertisements at the university.Before the tests were conducted, all subjects provided informed consent in accordance with the guidelines set by the institutional ethics review board.The study participants were organized into three groups of 20, sorted by BMI according to WHO standards (referenced in Table 1 of the manuscript) [10].The first group included individuals with a BMI under 25, identified as being within a healthy weight range.The second group, classified as overweight, consisted of participants with a BMI ranging from 25 to 30.The third group, defined as obese, included individuals with a BMI over 30.For the imaging process, all participants stood fully upright as X-ray images of their knee joints were taken in both the frontal and sagittal views using the DigitalDiagnost Rel 4.3 from Philips Medical Systems [11].A rheumatologist and an orthopedist thoroughly reviewed the X-ray images to verify that there was no apparent connection between osteoarthritis and obesity in any of the participants.

Data Collection
In this research, the external ground reaction forces and 3D movements of the lower limbs were measured (between 9 am and 1 pm) using the P6000 synchronized force platform and an optoelectronic motion capture system, both provided by BTS-Bioengineering, Inc.The motion capture equipment included eight SMART-DX EVO cameras, which were recording at a frequency of 100 Hz.We placed twenty-two spherical reflective markers, 20 mm in diameter, on vital anatomical landmarks, including the acromion, ASIS, sacrum, greater trochanter, femoral condyle, fibula head, lateral malleolus, fifth metatarsal head, and heel.Additional bar markers were affixed to the thigh and shank segments [12-14] .These markers, along with virtual markers established during periods of quiet standing, were employed to create anatomical coordinate systems for each lower limb segment [15, 16] .Participants walked barefoot at a self-chosen speed, completing at least five trials to ensure adequate data collection.The cutoff frequencies for the ground reaction forces and marker coordinates were 15 Hz and 6 Hz, respectively [17] .

Musculoskeletal model
An iterative kinematics-driven lower limb musculoskeletal model that accounts for the active-passive structures of the knee joint was developed.The hip and ankle were modeled as three-dimensional and two-dimensional spherical joints, respectively, surrounded by 31 muscles (4 around the ankle and 27 around the hip) [15,18].At the same time, the knee joint was

Joints surrounding muscles
Weight =77 kg) was scanned at Cleveland Clinic (Biomechanics laboratory) using a one Tesla extremity MRI scanner (Orthone, ONI Medical Systems Inc, Wilmington MA).A scanning protocol was used to provide a good contrast for the soft tissues within the same scan.The protocol characteristics are presented in Supplementary Table S(1).
Supplementary Table S(1): The magnetic resonance imaging settings (OpenKnee).The knee was placed in full extension, and the scanning process employed a 3D spoiled gradient-echo sequence with fat suppression, T.R. = 30, T.E.= 6.7,Flip Angle = 200, Field of View (FOV) = 150 mm X 150 mm, Slice Thickness = 1.5 mm.The imaging was conducted in three anatomical planes: (axial, sagittal, and coronal).About 18 minutes was spent finishing the process of scanning.These images are optimal to differentiate between the musculature, tendons, tissue fascia, and bone [19].The image data set was then imported into an MRI viewing and segmentation analysis package (3D slicer 4.8) and re-sampled in the sagittal, coronal, and axial planes.The muscle-bone junctions were identified from the MRI images following the procedure outlined in Dhaher and Kahn [20].Polygonal surfaces were used to generate a F.E. mesh of the knee joint using the Hypermesh (Altair Engineering, Troy, MI) and SOLIDWORKS (CAD) preprocessor.The structure of the tibiofemoral joint was adjusted to match the given dimension in the open knee public domain repository at Simtk.org [19].Bones were defined as rigid bodies [21] using 4-node quadrilateral elements that injunction with elastic boundaries with the articular cartilages.Eight-node hexahedral elements were used to represent the articular cartilages, (MPFL), quadriceps tendon (Q.T.), and patellar tendon (P.T.) cartilage layers and menisci are shown.More details on the system of axes and the joint center calculations can be found in [19].
The mesh of the model was obtained through a sensitivity analysis, where a maximum of 6% difference in the von-Mises stress was considered (STable.2).
Supplementary Table S(2): Mesh of the knee joint.

Fibril orientation
For each C3D8R element, a local coordinate system was used to define fibril orientation and implemented by a Python script.This script read the connectivity of each element and defined the local cross-sectional plane (x', y') and its normal vector along the local z'-direction (Sfig.4).

Supplementary Figure S(4):
The element local coordinate system (x', y', z') employed to define the orientation of the fibril.

Knee Model scaling
To match the representation of the subject's knee structure, the finite element (F.E.) knee model underwent morphing to align with the subject's anatomical measurements, as illustrated in Supplementary Figure S (5).This morphing process utilized anisotropic scaling, specifically by multiplying the nodal coordinates (x, y, z) by a percentage difference ratio.This ratio was calculated using the subject's measured dimensions, namely the maximum anterior-posterior

Interaction and loading analyses.
The explicit algorithm was used during all the simulations with a short step time (0.01s) to mimic quasi-static analysis.The boundary conditions were applied gradually using an exponential function for the normalized amplitude (Sfig.6) to attain smooth results.Frictionless interaction property was considered for driving surface-to-surface contact formulation.Computations were performed using an Intel(R) Core (T.M.) I9-12900KF@CPU3.20GHz,64.0 GB of RAM.

Supplementary Figure S(6):
The normalized amplitude as a function of the step time used to guide the application of boundary conditions during knee modeling.

Cartilage
The articular cartilage was modeled using incompressible hyperelastic fibrils reinforced composites behavior described by Sajjadinia et al., [23].The Cauchy stress (  ) in the used model was decomposed into a non-fibrillar (  ) and fibrillar (    ) parts as follow: Where F and J are the deformation gradient tensor and the volumetric deformation, respectively.n and   are the current direction and logarithmic strain of the fibril, respectively. 0 and   are the collagen stiffening coefficients (initial and strain-dependent).Gm is the shear modulus and  0  is an elastic material constant,   is the relative collagen fibril volume fraction.
The collagen networks were defined as primary and secondary bundles of fibrils based on their orientation relative to the articular cartilage depth (Sfig.7).The fibrils were oriented perpendicular to the subchondral junction and turned gradually in the middle zone to become parallel to the articular surface.For more details on the formulation of the material, please see prior works [23,24].A list of the properties of the material is presented in Supplementary Table S(3).Supplementary Table S(3): Articular cartilage materials properties.for the primary fibril and

Material parameters
for the secondary one. 2z: Normalized Depth of the articular cartilage (starting from the cartilage-bone junction area)

Meniscus
Since the isotropy of the transverse and axial plans in the meniscus has been thoroughly described, a particular subclass of orthotropy, transverse isotropy, was employed to represent the mechanical behavior of this substance [25][26][27].The local system axis of the meniscus is defined by axial, transverse, and circumferential axes, with the assumption that the transverse-axial plane is isotropic.As a result of this assumption, the number of independent constants in the matrix equals 5. To accomplish this, the transverse isotropy requires circumferential modulus (EC), transverse and axial modulus (Et=Ea), Poisson's ratio (νct= νca), which is defined as the ratio of the contractile strain in the transverse plane to the tensile strain in the circumferential direction under the load in the circumferential direction; Poisson's ratio νta, which is the Poisson's ratio within the transverse plane and shear modulus G (STable. 4

Ligaments
The knee tendons (PT and QT) were assumed to be neo-Hookean, with material coefficients (C10) of 55.9 MPa for the PT and 65.9 MPa for the QT [29].Meanwhile, knee ligaments were modeled using an incompressible transversely isotropic hyperelastic behavior [30] via an uncoupled representation of the strain energy function [31].The collagen fibers were uniformly distributed and properly bound to the isotropic and hyperelastic ground substance.The suggested strain energy function provides a silent reaction under any compressive loads and nonlinear stiffening behavior under tension as follows: Where   ,   and   are the strain energy's non-fibrilar, fibrillar, and volumetric parts.
n0 is the fiber orientation in the reference configuration, F is the deformation gradient tensor, c1, c2, c3 are the materials coefficients, and D is the incompressibility penalty parameter.The prestrains behavior was incorporated into the ligaments by decomposing the deformation gradient (F= F0 Fr) into a stress-free state (F0) and reference state (Fr).The pre-strains here were defined as the initial stretch ( 0 ) field with . The current study considered the sets of material parameters from our recent publication [28].
The predicted load-deflection curve of the tibial shaft exhibits a nonlinear stiffening behavior, as shown in Supplementary Figure S(9a).The computed axial deflection has a relatively higher increase at the initial equilibrium position under ligaments prestress and low axial load, then tend to stabilize at a higher load.On the side of the tibial plateaus load distribution, a higher load has been computed on the lateral compartment (~60%) than the medial one (~40%).The lateral plateau's covered area (meniscus cartilage) was the main supporter of the load, while on the medial plateau, almost the same proportions were computed between the covered and uncovered along the frontal axis while starting distally and switched at mid-range to proximally along the transversal axis (Sfig.10).These joint translations were computed based on the position of the primary tibial node (tibial bone reference node), which may alter if another node has been selected.
This may explain the slight computed differences with measured results under almost the same boundary conditions [37].In addition, the coupled internal rotation at the tibia substantially increased at high flexion angles.However, a small adduction rotation was computed during knee flexion and was characterized by an almost linear behavior.Finally, the overall kinematics of the tibiofemoral joint was in good agreement with the earlier cadaveric study of Germain et al., [37]. Supplementary

Muscles Optimization
A nonlinear optimization technique has been employed to evaluate the unknown muscle forces ({}) at each instance of the stance phase of gait.This optimization procedure minimizes an objective function of the sum of cubed muscle stresses (f (xi)) (4) under the constraint that muscle forces remain positive between their passive forces and total maximum active forces (6).All muscles' passive and maximum active forces were driven from a scaled musculoskeletal model that matched the subject's dimension considered in our gait analyses [15].The main constraint driving the evaluation of the muscle forces was the equilibrium equations (5).12).
represented by a complex nonlinear model consisting of full anatomical passive and active (8 muscles) structures (Sfig.1).Supplementary Figure S(1): Lower extremity musculoskeletal model, including the hip and ankle as spherical and hinge joints with an anterior and posterior view of the knee's three-dimensional finite element model.
An anatomically accurate model of knee joint consisting of the femur, tibial and patellar bones, articular surfaces, as well as the origins and insertions of the ligaments, was derived from digitized magnetic resonance image (MRI) transverse contours (OpenKnee public domain repository at Simtk.org).The knee specimen (female subject: Age 70 years; Height =170 cm; ligaments, and menisci (Sfig.2, 3).Supplementary Figure S(2): The workflow of building the 3D F.E. model of the knee joint; Clinical magnetic resonance image (MRI) of a knee joint, segmentation process, CAD file, and finite element model (F.E.) Supplementary Figure S(3): Anterior and posterior views of the three-dimensional finite element model of the knee showing the corresponding soft tissues and articular surfaces acting on the bones.Anterior cruciate ligaments (ACL), posterior cruciate ligaments (PCL), medial and lateral collateral ligament (MCL, LCL), lateral patellofemoral (LPFL), medial patellofemoral

(
MAP) length, maximum medial-lateral (MML) length, and the maximum tibiofemoral joint space width (MJW=max [MJW, LJW]), represented as (Dx, Dy, and Dz), in relation to the original model's dimensions [22].To execute the scaling, the mesh editing tool in Abaqus software was employed, coupled with a custom Python subroutine, enabling iterative adjustments to the model.Supplementary Figure S(5): measurements of the anatomical dimensions from the X-ray.
(7): Diagram showing the orientation and volume fraction of the collagen fibrils as a function of the depth in the articular cartilage.

(
cartilage-cartilage) areas.The total knee contact area increased from 404 to 1283 mm 2 when the tibial axial load was augmented by 2000 N (Sfig.9b).Average contact pressure on the tibial plateau reached its maximum value of 1.53 MPa under 2000 N axial load (Sfig.9c).Furthermore, the maximum compressive stress of 3.87 MPa, in the lateral plateau was computed in the uncovered zone near the central area.A stiffer moment-rotation response has been computed under increased pure coupled rotation (adduction and internal rotations).For example, the internal and adduction moments reach their maximum value of 5.6 and 29.3 Nm at 12 deg and 6 deg rotation, respectively (Sfig.9d).Our computed results under joint full extension posture have been corroborated by experimental measurements employing similar boundary conditions [33, 35, 38-41, 47-49].Supplementary Figure S(9): Tibial axial displacement (a), contact area (b), average and peak pressure (c, d) under axial compressive load up to 2000 N, tibial pure coupled momentrotation (e,f), along with experimental measurements under almost the same boundary conditions.The highlighted area indicated the standard deviation of experimental measurements.The tibial rolled on the femur posteriorly during passive knee flexion as the joint flexed to 90 degrees.The coupled joint laxities were directed laterally on the whole range of joint flexion

Figure S( 10 )
: Kinematics of the tibiofemoral joint during passive knee flexion.The highlighted area indicated the standard deviation of experimental measurements.(ac) tibial translation, (d-e) tibial coupled rotations.
For the forces in ligaments, The PCL force increased during knee flexion and reached its maximum value of 11 N at 90 degrees flexion.However, the ACL follows an opposite trend where a maximum of 77 N was computed at full extension and decreased after that on the rest of the range of joint flexion.The observed behavior confirms the primary role of the cruciate ligaments, in which the ACL is the primary component that resists anterior tibial translation at full extension.At the same time, the PCL remained slack and got active at higher angle flexion to resist associated posterior tibial translation.The following figure (Sfig.11) compares computed cruciate ligament forces with computational prediction and experimental measurements under almost the same boundary conditions [43-45, 50].Supplementary Figure S(11): Comparison of our predicted and measured cruciate ligament forces during passive knee flexion.The highlighted area indicated the standard deviation of experimental measurements.
{  } ≤ {} ≤ {  } (6) With   ,   , [R] ,   ,   are the force of muscle, a passive force of muscle, and lever arms matrix at different instances during the stance phase, maximum force, and physiological cross-sectional areas for muscle i, respectively [15].[] are the (required) lower limb joint moments (hip, knee, and ankle) computed during the different instances of the stance phase (Sfig.

(femur, tibia and Patella) 27560 S4R
A list of the properties of the material is presented in supplementary Table S(4).